In [24]:

    

import numpy 
import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt


%matplotlib inline
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16


    

In [25]:

    

nx=100                           # the number of nodes in lattice direction 
dt=1
dx=1
x1=0.
x2=1.0                           # the length 
alpha=0.25                       # heat diffusion coefficient 
D=1                              # the dimension of the problem 
mstep=6000                       # the number of time step
omega=1./(0.5+(alpha*D))         #Chapman-Enskog  dt=1 and dx=1 
Tleft=1.0                        # left wall temperature
Tright=0.0                       # right wall temperature
k=3                              # k=0,1,2    <==== c(2)===c(0)====> c(1)


    

In [26]:

    

x=numpy.linspace(x1,x2,nx+1)
w=numpy.zeros(k)                           # weitghting factor
T=numpy.zeros(nx+1)                       # Temperature matrix
f= numpy.zeros((k, nx+1))                 # distribution function
feq= numpy.zeros((k, nx+1))               # Equilibrum distribution function


    

In [27]:

    

w[0]=0.0
w[1]=0.5
w[2]=0.5  # k=3


    

In [28]:

    

T[:]=0.0                                   #initial temperature
T[0]=1.0
for i in range (k): #k=0,1,2
 f[i,:]=w[i]*T[:]


    

In [29]:

    

#================== Initial value 

T[:]=0.0 #initial temperature
T[0]=1.0
for i in range (k): #k=0,1,2
 f[i,:]=w[i]*T[:]


    

In [30]:

    

#   Main loop  : comprised two parts :collision and streaming

#=====================
for n in range(mstep) :
                                                        # collision process
 for i in range(nx+1):
  T[i]=f[0,i]+f[1,i]+f[2,i]
  feq[0:k,i]=w[0:k]*T[i]
  f[0:k,i]=(1.-omega)*f[0:k,i]+omega*feq[0:k,i]
                                                        # streaming process
 for i in range(0,nx):
  f[1,nx-i]=f[1,nx-i-1]
  f[2,i]=f[2,i+1]

                                                    # Boundary condition 
                        
 f[1,0]=Tleft-f[2,0]   # constant temperature at left  T=1,0
 f[2,nx]=Tright-f[1,nx] # constant temperature at right T=0.0
 T[0]=1.0
 T[nx]=0.0
                        
# end of the main loop


    

In [31]:

    

#     Finite difference # 


Tf=numpy.zeros(nx+1)   # Temperature finite difference
To=numpy.zeros(nx+1)   # Temperature storaage for previous time
Tf[0]=1.0
Tf[nx]=0.0
dx=1./nx
dt=1e-4                    #  2*alpha*dt/dx^2  << 1 === > for stability
for n in range(mstep) :
 To[:]=Tf[:]
 for i in range(1,nx):
  Tf[i]=To[i]+ ( (dt*alpha/(dx**2.)) *(To[i+1]-2*To[i]+To[i-1]) )


    

In [34]:

    

plt.figure(figsize=(15,6), dpi=100)
plt.xlabel('x', fontsize=10) #x label
plt.ylabel('T', fontsize=10) #y label
plt.plot(x,T, 'r-',label=' Lattice Boltzmann Method'); 
plt.plot(x,Tf, 'bo',label=' Finite Difference Method '); 
plt.legend();


    

In [ ]: